Lab 11

Author

Erin Yu

Lab 11

Read in Data

library(data.table)


cv_states <- as.data.frame(fread("https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv"))


state_pops <- as.data.frame(fread("https://raw.githubusercontent.com/COVID19Tracking/associated-data/master/us_census_data/us_census_2018_population_estimates_states.csv"))


state_pops$abb <- state_pops$state
state_pops$state <- state_pops$state_name
state_pops$state_name <- NULL


cv_states <- merge(cv_states, state_pops, by = "state")

eval(TRUE)
[1] TRUE

Look at Data

dim(cv_states)
[1] 58094     9
head(cv_states)
    state       date fips   cases deaths geo_id population pop_density abb
1 Alabama 2023-01-04    1 1587224  21263      1    4887871    96.50939  AL
2 Alabama 2020-04-25    1    6213    213      1    4887871    96.50939  AL
3 Alabama 2023-02-26    1 1638348  21400      1    4887871    96.50939  AL
4 Alabama 2022-12-03    1 1549285  21129      1    4887871    96.50939  AL
5 Alabama 2020-05-06    1    8691    343      1    4887871    96.50939  AL
6 Alabama 2021-04-21    1  524367  10807      1    4887871    96.50939  AL
tail(cv_states)
        state       date fips  cases deaths geo_id population pop_density abb
58089 Wyoming 2022-09-11   56 175290   1884     56     577737    5.950611  WY
58090 Wyoming 2022-08-21   56 173487   1871     56     577737    5.950611  WY
58091 Wyoming 2021-01-26   56  51152    596     56     577737    5.950611  WY
58092 Wyoming 2021-02-21   56  53795    662     56     577737    5.950611  WY
58093 Wyoming 2021-08-22   56  70671    809     56     577737    5.950611  WY
58094 Wyoming 2021-03-20   56  55581    693     56     577737    5.950611  WY
str(cv_states)
'data.frame':   58094 obs. of  9 variables:
 $ state      : chr  "Alabama" "Alabama" "Alabama" "Alabama" ...
 $ date       : IDate, format: "2023-01-04" "2020-04-25" ...
 $ fips       : int  1 1 1 1 1 1 1 1 1 1 ...
 $ cases      : int  1587224 6213 1638348 1549285 8691 524367 1321892 1088370 1153149 814025 ...
 $ deaths     : int  21263 213 21400 21129 343 10807 19676 16756 16826 15179 ...
 $ geo_id     : int  1 1 1 1 1 1 1 1 1 1 ...
 $ population : int  4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 ...
 $ pop_density: num  96.5 96.5 96.5 96.5 96.5 ...
 $ abb        : chr  "AL" "AL" "AL" "AL" ...

There are 9 variables including “state,” “date,” “fips,” “cases,” “deaths,” “geo_id,” “population,” “pop_density,” and “abb” with a total of 58,094 observations. The first 5 observations are from Alabama with number of deaths ranging from 213 to 21,400. The last 5 observations are from Wyoming, with number of deaths ranging from 596 to 1,884.

Format the data

cv_states$date <- as.Date(cv_states$date, format="%Y-%m-%d")


state_list <- unique(cv_states$state)
cv_states$state <- factor(cv_states$state, levels = state_list)
abb_list <- unique(cv_states$abb)
cv_states$abb <- factor(cv_states$abb, levels = abb_list)


cv_states <- cv_states[order(cv_states$state, cv_states$date),]


str(cv_states)
'data.frame':   58094 obs. of  9 variables:
 $ state      : Factor w/ 52 levels "Alabama","Alaska",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ date       : Date, format: "2020-03-13" "2020-03-14" ...
 $ fips       : int  1 1 1 1 1 1 1 1 1 1 ...
 $ cases      : int  6 12 23 29 39 51 78 106 131 157 ...
 $ deaths     : int  0 0 0 0 0 0 0 0 0 0 ...
 $ geo_id     : int  1 1 1 1 1 1 1 1 1 1 ...
 $ population : int  4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 4887871 ...
 $ pop_density: num  96.5 96.5 96.5 96.5 96.5 ...
 $ abb        : Factor w/ 52 levels "AL","AK","AZ",..: 1 1 1 1 1 1 1 1 1 1 ...
head(cv_states)
       state       date fips cases deaths geo_id population pop_density abb
1029 Alabama 2020-03-13    1     6      0      1    4887871    96.50939  AL
597  Alabama 2020-03-14    1    12      0      1    4887871    96.50939  AL
282  Alabama 2020-03-15    1    23      0      1    4887871    96.50939  AL
12   Alabama 2020-03-16    1    29      0      1    4887871    96.50939  AL
266  Alabama 2020-03-17    1    39      0      1    4887871    96.50939  AL
78   Alabama 2020-03-18    1    51      0      1    4887871    96.50939  AL
tail(cv_states)
        state       date fips  cases deaths geo_id population pop_density abb
57902 Wyoming 2023-03-18   56 185640   2009     56     577737    5.950611  WY
57916 Wyoming 2023-03-19   56 185640   2009     56     577737    5.950611  WY
57647 Wyoming 2023-03-20   56 185640   2009     56     577737    5.950611  WY
57867 Wyoming 2023-03-21   56 185800   2014     56     577737    5.950611  WY
58057 Wyoming 2023-03-22   56 185800   2014     56     577737    5.950611  WY
57812 Wyoming 2023-03-23   56 185800   2014     56     577737    5.950611  WY
head(cv_states)
       state       date fips cases deaths geo_id population pop_density abb
1029 Alabama 2020-03-13    1     6      0      1    4887871    96.50939  AL
597  Alabama 2020-03-14    1    12      0      1    4887871    96.50939  AL
282  Alabama 2020-03-15    1    23      0      1    4887871    96.50939  AL
12   Alabama 2020-03-16    1    29      0      1    4887871    96.50939  AL
266  Alabama 2020-03-17    1    39      0      1    4887871    96.50939  AL
78   Alabama 2020-03-18    1    51      0      1    4887871    96.50939  AL
summary(cv_states)
           state            date                 fips           cases         
 Washington   : 1158   Min.   :2020-01-21   Min.   : 1.00   Min.   :       1  
 Illinois     : 1155   1st Qu.:2020-12-06   1st Qu.:16.00   1st Qu.:  112125  
 California   : 1154   Median :2021-09-11   Median :29.00   Median :  418120  
 Arizona      : 1153   Mean   :2021-09-10   Mean   :29.78   Mean   :  947941  
 Massachusetts: 1147   3rd Qu.:2022-06-17   3rd Qu.:44.00   3rd Qu.: 1134318  
 Wisconsin    : 1143   Max.   :2023-03-23   Max.   :72.00   Max.   :12169158  
 (Other)      :51184                                                          
     deaths           geo_id        population        pop_density       
 Min.   :     0   Min.   : 1.00   Min.   :  577737   Min.   :    1.292  
 1st Qu.:  1598   1st Qu.:16.00   1st Qu.: 1805832   1st Qu.:   43.659  
 Median :  5901   Median :29.00   Median : 4468402   Median :  107.860  
 Mean   : 12553   Mean   :29.78   Mean   : 6397965   Mean   :  423.031  
 3rd Qu.: 15952   3rd Qu.:44.00   3rd Qu.: 7535591   3rd Qu.:  229.511  
 Max.   :104277   Max.   :72.00   Max.   :39557045   Max.   :11490.120  
                                                     NA's   :1106       
      abb       
 WA     : 1158  
 IL     : 1155  
 CA     : 1154  
 AZ     : 1153  
 MA     : 1147  
 WI     : 1143  
 (Other):51184  
min(cv_states$date)
[1] "2020-01-21"
max(cv_states$date)
[1] "2023-03-23"

The date range is from 2020-01-21 to 2023-03-23. The range of cases are from 1 to 12,169,158. The range of deaths is 0 to 104,277.

Add new_cases and new_deaths and correct outliers

library(ggplot2)
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:data.table':

    between, first, last
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(plotly)

Attaching package: 'plotly'
The following object is masked from 'package:ggplot2':

    last_plot
The following object is masked from 'package:stats':

    filter
The following object is masked from 'package:graphics':

    layout
library(zoo)

Attaching package: 'zoo'
The following objects are masked from 'package:data.table':

    yearmon, yearqtr
The following objects are masked from 'package:base':

    as.Date, as.Date.numeric
for (i in 1:length(state_list)) {
  cv_subset = subset(cv_states, state == state_list[i])
  cv_subset = cv_subset[order(cv_subset$date),]

 
  cv_subset$new_cases = 0
  cv_subset$new_deaths = 0

  for (j in 2:nrow(cv_subset)) {
    cv_subset$new_cases[j] = cv_subset$cases[j] - cv_subset$cases[j-1]
    cv_subset$new_deaths[j] = cv_subset$deaths[j] - cv_subset$deaths[j-1]
  }

 
  cv_states$new_cases[cv_states$state==state_list[i]] = cv_subset$new_cases
  cv_states$new_deaths[cv_states$state==state_list[i]] = cv_subset$new_deaths
}


cv_states <- cv_states %>% dplyr::filter(date >= "2021-06-01")


p1 <- ggplot(cv_states, aes(x = date, y = new_cases, color = state)) + geom_line() + geom_point(size = 0.5, alpha = 0.5)
ggplotly(p1)
p1 <- NULL 

p2 <- ggplot(cv_states, aes(x = date, y = new_deaths, color = state)) + geom_line() + geom_point(size = 0.5, alpha = 0.5)
ggplotly(p2)
p2 <- NULL 


cv_states$new_cases[cv_states$new_cases < 0] = 0
cv_states$new_deaths[cv_states$new_deaths < 0] = 0


for (i in 1:length(state_list)) {
  cv_subset = subset(cv_states, state == state_list[i])

  cv_subset$cases[1] = cv_subset$cases[1]
  cv_subset$deaths[1] = cv_subset$deaths[1]

  for (j in 2:nrow(cv_subset)) {
    cv_subset$cases[j] = cv_subset$new_cases[j] + cv_subset$cases[j-1]
    cv_subset$deaths[j] = cv_subset$new_deaths[j] + cv_subset$deaths[j-1]
  }


  cv_states$cases[cv_states$state == state_list[i]] = cv_subset$cases
  cv_states$deaths[cv_states$state == state_list[i]] = cv_subset$deaths
}


cv_states$new_cases = zoo::rollmean(cv_states$new_cases, k=7, fill=NA, align='right') %>% round(digits = 0)
cv_states$new_deaths = zoo::rollmean(cv_states$new_deaths, k=7, fill=NA, align='right') %>% round(digits = 0)


p2 <- ggplot(cv_states, aes(x = date, y = new_deaths, color = state)) + geom_line() + geom_point(size = 0.5, alpha = 0.5)
ggplotly(p2)
p2 <- NULL

Add additional variables

cv_states$per100k = as.numeric(format(round(cv_states$cases / (cv_states$population / 100000), 1), nsmall = 1))
cv_states$newper100k = as.numeric(format(round(cv_states$new_cases / (cv_states$population / 100000), 1), nsmall = 1))
Warning: NAs introduced by coercion
cv_states$deathsper100k = as.numeric(format(round(cv_states$deaths / (cv_states$population / 100000), 1), nsmall = 1))
cv_states$newdeathsper100k = as.numeric(format(round(cv_states$new_deaths / (cv_states$population / 100000), 1), nsmall = 1))
Warning: NAs introduced by coercion
cv_states = cv_states %>% mutate(naive_CFR = round((deaths * 100 / cases), 2))


cv_states_today = subset(cv_states, date == max(cv_states$date))

Scatterplots

Create a scatterplot usng plot_ly

library(plotly)


plot_cases <- cv_states_today %>% 
  plot_ly(x = ~population, y = ~cases, 
          type = 'scatter', mode = 'markers', color = ~state,
          size = ~population, sizes = c(5, 70), marker = list(sizemode='diameter', opacity=0.5)) %>%
  layout(title = "Population vs. Cases",
         xaxis = list(title = "Population"),
         yaxis = list(title = "Cases"))
plot_cases
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
cv_states_today_filter <- cv_states_today %>% filter(state != "District of Columbia")


plot_cases_filtered <- cv_states_today_filter %>% 
  plot_ly(x = ~population, y = ~cases, 
          type = 'scatter', mode = 'markers', color = ~state,
          size = ~population, sizes = c(5, 70), marker = list(sizemode='diameter', opacity=0.5)) %>%
  layout(title = "Population vs. Cases (Filtered)",
         xaxis = list(title = "Population"),
         yaxis = list(title = "Cases"))
plot_cases_filtered
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
plot_deaths_per100k <- cv_states_today_filter %>% 
  plot_ly(x = ~population, y = ~deathsper100k, 
          type = 'scatter', mode = 'markers', color = ~state,
          size = ~population, sizes = c(5, 70), marker = list(sizemode='diameter', opacity=0.5)) %>%
  layout(title = "Population vs. Deaths per 100k",
         xaxis = list(title = "Population"),
         yaxis = list(title = "Deaths per 100k"))
plot_deaths_per100k
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
plot_hover <- cv_states_today_filter %>% 
  plot_ly(x = ~population, y = ~deathsper100k, 
          type = 'scatter', mode = 'markers', color = ~state,
          size = ~population, sizes = c(5, 70), marker = list(sizemode='diameter', opacity=0.5),
          hoverinfo = 'text',
          text = ~paste(state, "<br>Cases per 100k:", per100k, "<br>Deaths per 100k:", deathsper100k)) %>%
  layout(title = "Population-normalized COVID-19 deaths (per 100k) vs. Population",
         xaxis = list(title = "Population"),
         yaxis = list(title = "Deaths per 100k"),
         hovermode = "compare")
plot_hover
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors

Explore scatterplot trend interactively using ggplotly and geom_smooth

library(ggplot2)
library(plotly)

cv_states_today_scatter <- cv_states_today_filter
p <- ggplot(cv_states_today_scatter, aes(x = population, y = newdeathsper100k, size = population)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", color = "blue") +
  labs(title = "Population Density vs New Deaths per 100k",
       x = "Population Density",
       y = "New Deaths per 100k")

ggplotly(p)
`geom_smooth()` using formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: size.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
  the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
  variable into a factor?

The number of new deaths per 100k seem to decrease with increasing population density. It seems that there may be a negative correlation.

Multiple Line Chart

plot_naive_cfr <- plot_ly(cv_states, x = ~date, y = ~naive_CFR, color = ~state, type = "scatter", mode = "lines") %>%
  layout(title = "Naive Case Fatality Rate (CFR) Over Time",
         xaxis = list(title = "Date"),
         yaxis = list(title = "Naive CFR (%)"))
plot_naive_cfr
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
Returning the palette you asked for with that many colors
plot_florida <- cv_states %>% filter(state == "Florida") %>%
  plot_ly(x = ~date, y = ~new_cases, name = "New Cases", type = "scatter", mode = "lines", line = list(color = "blue")) %>%
  add_trace(y = ~new_deaths, name = "New Deaths", mode = "lines", line = list(color = "red")) %>%
  layout(title = "New Cases and New Deaths in Florida",
         xaxis = list(title = "Date"),
         yaxis = list(title = "Count"),
         hovermode = "x unified")
plot_florida

The states that increased in september included South Carolina, Florida, and Oklahoma. Over the next 2 months they seemed to gradually increase. in naive_cfr. The time difference between the peak cases and deaths in Florida was about half a month

Heat maps

library(tidyr)
cv_states_mat <- cv_states %>% select(state, date, new_cases) %>% dplyr::filter(date > as.Date("2021-06-15"))
cv_states_mat2 <- as.data.frame(pivot_wider(cv_states_mat, names_from = state, values_from = new_cases))
rownames(cv_states_mat2) <- cv_states_mat2$date
cv_states_mat2$date <- NULL
cv_states_mat2 <- as.matrix(cv_states_mat2)

plot_ly(x = colnames(cv_states_mat2), y = rownames(cv_states_mat2),
        z = ~cv_states_mat2,
        type = "heatmap",
        showscale = T)
cv_states_mat <- cv_states %>% select(state, date, newper100k) %>% dplyr::filter(date > as.Date("2021-06-15"))
cv_states_mat2 <- as.data.frame(pivot_wider(cv_states_mat, names_from = state, values_from = newper100k))
rownames(cv_states_mat2) <- cv_states_mat2$date
cv_states_mat2$date <- NULL
cv_states_mat2 <- as.matrix(cv_states_mat2)

plot_ly(x = colnames(cv_states_mat2), y = rownames(cv_states_mat2),
        z = ~cv_states_mat2,
        type = "heatmap",
        showscale = T)
filter_dates <- seq(as.Date("2021-06-15"), max(cv_states$date), by = "2 weeks")
cv_states_mat <- cv_states %>% select(state, date, newper100k) %>% filter(date %in% filter_dates)
cv_states_mat2 <- as.data.frame(pivot_wider(cv_states_mat, names_from = state, values_from = newper100k))
rownames(cv_states_mat2) <- cv_states_mat2$date
cv_states_mat2$date <- NULL
cv_states_mat2 <- as.matrix(cv_states_mat2)

plot_ly(x = colnames(cv_states_mat2), y = rownames(cv_states_mat2),
        z = ~cv_states_mat2,
        type = "heatmap",
        showscale = T)

California, New York, Texas, and Washington DC stood out from the heatmap of new cases. When repeated with newper100k, Rhode Island stood out the most.

Map

library(dplyr)
library(plotly)


pick.date <- "2021-10-15"


cv_per100_date <- cv_states %>% 
  filter(date == pick.date) %>% 
  select(state, abb, newper100k, cases, deaths) %>%
  mutate(
    state_name = state,
    state = abb,
    hover = paste(state_name, '<br>', 
                  "Cases per 100k: ", newper100k, '<br>', 
                  "Cases: ", cases, '<br>', 
                  "Deaths: ", deaths)
  )


cv_per100_today <- cv_states_today %>% 
  select(state, abb, newper100k, cases, deaths) %>%
  mutate(
    state_name = state,
    state = abb,
    hover = paste(state_name, '<br>', 
                  "Cases per 100k: ", newper100k, '<br>', 
                  "Cases: ", cases, '<br>', 
                  "Deaths: ", deaths)
  )


set_map_details <- list(
  scope = 'usa',
  projection = list(type = 'albers usa'),
  showlakes = TRUE,
  lakecolor = toRGB('white')
)


shadeLimit <- 125


fig_pick_date <- plot_geo(cv_per100_date, locationmode = 'USA-states') %>% 
  add_trace(
    z = ~newper100k, text = ~hover, locations = ~state,
    color = ~newper100k, colors = 'Purples'
  ) %>%
  colorbar(title = paste0("Cases per 100k: ", pick.date), limits = c(0, shadeLimit)) %>%
  layout(
    title = paste('Cases per 100k by State as of ', pick.date, '<br>(Hover for value)'),
    geo = set_map_details
  )


fig_today <- plot_geo(cv_per100_today, locationmode = 'USA-states') %>% 
  add_trace(
    z = ~newper100k, text = ~hover, locations = ~state,
    color = ~newper100k, colors = 'Purples'
  ) %>%
  colorbar(title = paste0("Cases per 100k: ", Sys.Date()), limits = c(0, shadeLimit)) %>%
  layout(
    title = paste('Cases per 100k by State as of ', Sys.Date(), '<br>(Hover for value)'),
    geo = set_map_details
  )


subplot(fig_pick_date, fig_today, nrows = 2, margin = 0.05)

The CFR has changed from some medium to dark purple spots on the map from 10/15/2021 to nearly all white on today’s date